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Abstract Structure and ordering in swarms of active particles have much in com¬ 
mon with condensed matter systems like magnets or liquid crystals. A number 
of important characteristics of such materials can be obtained via dynamic tests 
such as hysteresis. In this work, we show that dynamic hysteresis can be observed 
also in swarms of active particles and possesses similar properties to the coun¬ 
terparts in magnetic materials. To study the swarm dynamics, we use computer 
simulations of the active Brownian particle model with dissipative interactions. 
The swarm is confined to a narrow linear channel and the one-dimensional polar 
order parameter is measured. In an oscillating external field, the order parameter 
demonstrates dynamic hysteresis with the shape of the loop and its area varying 
with the amplitude and frequency of the applied field, swarm density and the noise 
intensity. We measure the scaling exponents for the hysteresis loop area, which can 
be associated with the controllability of the swarm. Although the exponents are 
non-universal and depend on the system’s parameters, their limiting values can 
be predicted using a generic model of dynamic hysteresis. We also discuss similar¬ 
ities and differences between the swarm ordering dynamics and two-dimensional 
magnets. 

PACS 05.65.+b • 64.70.qj ■ 87.18.Nq 


1 Introduction 

Hysteresis is a nonlinear phenomenon commonly observed in various metastable 
systems, which have more than one internal state [I:. In such systems, the response 
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to a change of the environment depends on the history of this environment. Al¬ 
though the most prominent example of hysteresis is a magnetisation response of 
a ferromagnet in an oscillating field [2], this effect is not limited to only ferro¬ 
magnetic or ferroelectric materials. Hysteresis is also observed in various physical, 
mechanical, chemical, biological and ecological systems. In active systems [3MII5] 
hysteresis of the collective motion states has been shown to arise near the point 
of orientational phase transition, where system is most sensitive to any changes in 
fluctuation strength. 

While the interest in hysteretic effects is driven mostly by practical applications 
in electronics and engineering, hysteresis occurring in many-body systems is an 
intriguing fundamental problem on itself [Bim . Although the phenomenon has been 
known for several hundreds of years, the systematic study of hysteresis has started 
only in the last quarter of the twentieth century [T|. 

In this study, we explore the dynamic hysteresis in swarms of active species 
using the magnetic analogy. This analogy is based on the ability of active swarms 
to reach orientationally ordered states, similar to those in ferromagnets |8H9llT0l im 
mmm- The transition from a disordered state of a swarm of interacting active 
particles to an ordered state happens upon reduction of noise at fixed propulsion 
speed. The common polar order parameter for the swarm, the mean particle veloc¬ 
ity, measured as a function of the noise amplitude, behaves in the same way as the 
magnetisation vector in magnetic materials upon the variation of temperature m 
nn. The similarity between swarms and magnets is not limited to the behaviour of 
the mean order parameter but covers a wide range of more subtle properties such 
as the spatial correlations (two-point correlation function), susceptibility, and the 
divergence of the correlation radius at the critical point mmmm- The ideas 
from lattice models of magnets such as Heisenberg model and Ising model have 
been successfully applied to describe swarm dynamics |10U18| . 

We hope to extend the analogy between the equilibrium condensed matter and 
active swarms to dynamic properties. Although the interaction of active swarms 
with external fields has not been extensively studied so far, it is easy to envision 
that an orienting field would have the same effect on a swarm that is observed in 
magnetic systems. In particular, as the swarm has a finite orientation relaxation 
time, there must be a room for dynamic hysteresis controlled by a competition 
of the external drive and the internal relaxation. Therefore, one can attempt to 
characterise the swarm dynamics by retentivity , the ability to align in absence of 
the field, coercivity , the magnitude of field in the opposite direction needed to 
revert the direction of motion of the swarm, and susceptibility , the intensity of the 
response of the swarm to the action of field. In addition, one can hope to learn 
the main relaxation times, for example the time needed to reorient the swarm, 
and estimate the strength and the frequency of the required controlling signals. 
All these quantities can become extremely useful if we try to control the collective 
dynamics in either robotic swarms, human crowds or groups of animals in a farm 
or in Nature. 

The remainder of the paper is organised as follows: Section [2] describes the 
active Brownian particle model with interactions and simulation settings, Section 
[3] presents the results on orientational ordering in the active swarm, Section [4] 
presents the discussion of our main observations, and Section [5] concludes the 
paper. 





Orientational hysteresis in active swarms 


3 


2 Model and simulation setup 


To study the orientational hysteresis in an active swarm, we use a two-dimensional 
system of active Brownian particles (ABP) |19ll20| with dissipative interactions. 
Motion of the agents is confined into a narrow straight channel with periodic 
boundaries in the long direction and purely repulsive walls in the short direction. 
The ABP-DI model is able to produce a globally aligned phase if the energy influx 
rate is sufficiently high and the interactions are sufficiently strong ED- 

Our implementation of the ABP-DI is as described in our previous works EH 
f22j . The particle motion is governed by the Langevin equation 



( 1 ) 


where m is the particle mass (set to unity in this work) and V,; is the velocity of 
particle i. The total force F i(t) acting on each particle is given by: 


F, =Ff - 7 B V, + Ff + V2D%(f) + H(i) (2) 


where Ff is the force that comes from interactions within the swarm, 7 E is the 
coefficient of viscous friction, which is set by the properties of the environment, Ff 
is the thrust term. The term %/2 D E £i(t) is the random force of strength D E and 
£(i) is representing a Gaussian white noise with zero-mean and unit variance. The 
strength of the noise is set by the fluctuation-dissipation relation at the ambient 
temperature T E 
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The thrust term has the form: 


F 7 = 


dq 


c + dV/ 


rV, 


( 4 ) 


where d is the constant determining the rate of conversion of internal energy into 
kinetic energy, c is the parameter controlling energy loss, and q is the constant 
determining the gain of energy from the environment. In a stationary state, mo¬ 
tion of each particle is characterized by velocity P 0 2 , which is defined through the 
system’s parameters as 


Vo 


Q c 
' y E d 


( 5 ) 


at q > ^ E c/d p~9l[23] . 

Ff (t) is the force coming from the interactions within the swarm, which de¬ 
scribes inelastic collisions between the particles according to the Dissipative Par¬ 
ticle Dynamics (DPD) method EH- 

The dissipative force is taken in the form of a friction force applied to the 
component of the motion in the direction of the particle connecting vector. It 
generally consists of three parts: 


£(Fg + Fg + Fg) 


(6) 
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where F^, F-, and F^ represent the conservative, dissipative, and random forces 
between particles i and j, respectively. The conservative force reflects the excluded 
volume interactions: 



Tij < r r 
r%j > r r 


(7) 


where r. (;/ = rj — r j is the distance between particles i and j, r t j = |r ^ | is its magni¬ 
tude, Tij = rij/rij is the unit vector from j to i, G is a parameter determining the 
maximum repulsion between the particles, r r = r c /2 is the radius of the repulsion 
zone, and r c is the cut-off distance of the interaction. 

The dissipative force punishes the velocity differences between the neighbouring 
particles and, therefore, provides a mechanism of relaxation of the velocity field 
towards the stationary state. We take it in the form of a friction force applied to 
the component of the relative motion in the direction of the particle connecting 
vector ED, i.e. a velocity adjustment for particles following one another. 

F ?j = -7 S u(r ij )(r ij ■ Vijjfy (8) 

The pairwise friction coefficient 7 s determines the degree of inelasticity of the 
collisions, and co(r) is a weight function that describes the particle’s “soft shell”: 

cH.K'-y-) 2 ' (9) 

( 0 , rij > r c 

We neglect the random pairwise force Ffj in this work. 

The external field H(f) in our model is set by 


F[(t) = Hg sin ( u>t)x 


( 10 ) 


where Ho is the amplitude of the periodic field, to = 2nf is the angular frequency, 
and x is the unit vector pointing along the s-axis. 

Since the periodic force acts only along one axis, the natural order parameter 
for our system is the mean one-component velocity V x 

1 N 

= ( n ) 

i =1 


We used non-normalised order parameter in the study of the hysteresis loops. To 
study the phase behaviour of the swarm in absence of the field, we normalised the 
order parameter to bring it into the range from 0 to 1 as follows 


(fin — 


\ Vi(t) /’ 


( 12 ) 


where (•) stands for ensemble average. To locate the phase transition points pre¬ 
cisely we calculated the Binder cumulant ES using the normalised order parameter 
defined by Eq. (1121) 

(P L )t 

3K>? 


g l = 1 


(13) 
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(a) (b) 



Fig. 1 Dynamics of the orientational order parameter ip in the ABP-DI system subjected to 
a constant field (p = 0.04, T E = 0.3). The field is switched on at to = 800. (a) Varying Ho, 
q = 1 . (b) Varying q, Hq = 0.2. 


where (-)t stands for the time average and index L denotes the value calculated 
for a system of linear length L. The Binder cumulant has a very weak dependence 
on the system size so Gl takes a universal value at the critical point for any L if 
the density p and the energy influx rate q are kept constant [26]. In this work, the 
transition points in q — p plane are, therefore, defined as an intersection between 
the three Gl — q curves for different channel length L at constant density. Those 
points were then used to construct the phase diagram. 

The area of hysteresis loops has been calculated as 


A = 


<p(H)dH 


(14) 


The system relaxation time r was measured from the order parameter relaxation 
dynamics towards the steady state ipoo at fixed p, q, T E , and Hq 


<p(t) = (poo (l - e (t to)/T ). 
upon application of a step-like signal at time £q: 


H(t) 


Hq, t > to 
0 , t < to 


(15) 


(16) 


Examples of the order parameter relaxation in constant field are shown in Fig. 
U| Dependence of the relaxation time r on the strength of applied field Ho and 
temperature T E is presented in Fig. [2] 

All simulations were performed with the following set of key parameters: = 

0.3, G = 1, d = 2, c = 0.8, 7 s = 0.3, r c = 2 (cut-off radius), r T = 1 (the particle 
radius of repulsion). The radius of repulsion of the particles sets the unit of length 
if the simulation system. To set the unit of time, we choose a unit speed v = 1 
such that a particle moving at V = v would make a distance r r per unit time. This 
definition can be reformulated in terms of kinetic energy: our simulation units are 
such that an active particle moving at a speed of one body radius per unit time 
would have a kinetic energy E = mV 2 /2 = 1/2. Therefore, a temperature T E = 0.3 
in our settings, which sets the noise amplitude, means that the root-mean-square 
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(a) (b) 




Fig. 2 Relaxation time of the orientational order parameter tp in the ABP-DI system subjected 
to a constant field (p = 0.04, q = 1). Each point represents an average over five independent 
runs, (a) Varying Hq , T E = 0.3 . (b) Varying T E , Hq = 1. 


speed of particles without propulsion (q = 0) is Vrms = \/T E /m = 0.548, i.e. 0.548 
body radii per unit time. 

Simulations were performed with time step of At = 0.01. The positions of the 
agents were propagated using the Verlet algorithm m- The geometric confinement 
was represented by the linear channel with its walls lying along the a:-axis and 
periodic boundaries in the z-direction. Dimensions of the channel were fixed at 50 x 
500 units for all runs (unless stated otherwise). Repulsions from the channel walls 
were modeled as mirror-like reflections: after bouncing off the wall the component 
of the particle velocity normal to the wall was getting the opposite sign while 
the tangential one was kept unchanged. In all our simulations the oscillating field 
was applied parallel to the confining walls. Initial positions of particles have been 
chosen at random in all simulations. Total number of time steps has been set in 
the region 2 x 10® -lx 10 8 depending on the frequency of the external field and 
the hysteresis loops were averaged over at least 10 cycles. 


3 Results 

3.1 Orientational ordering 

As we found in our previous work, an ABP-DI swarm is known to order orienta- 
tionally in free space at sufficiently high density, strong interactions (as expressed 
by F s ), and/or strong propulsion [21]. We observe the same behaviour in confine¬ 
ment. A few typical snapshots of the part of the system are shown in Figs. [3] and 
[4] The distribution of the active particles along and across the channel is visibly 
affected by their incoming power q (see Fig. [3]) and by the number density p (see 
Fig. [4]). While at low incoming energy rates the particles behave like a gas with not 
much correlation in their motion, at certain critical level of energy pumping their 
motion becomes orientationally ordered. Another obvious result of increasing the 
incoming power is the particle aggregation. As we see in Fig. (3 }d, already at q = 1 
we observe significant density fluctuations and at q = 10 large compact clusters 
appear. A similar effect is observed on increasing particle number density p: at low 
density p = 0.04 (Fig. [Ur) we see only gas-like behaviour and disordered motion, 
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Fig. 3 Simulation snapshots (cutout of the full simulation box) of swarms for the ABP-DI 
model in a linear channel confinement at various energy input rate q (p = 0.1, T E = 0.3). (a) 
q = 0. (b) q = 1. (c) q = 10. Arrows indicate the direction of motion of the individuals as well 
as the magnitudes of the velocities. 


while at higher densities p = 0.1 and p = 0.4 both the aggregation and alignment 
become pronounced. 

In Fig. [5^, we plot the polar order parameter as a function of the energy influx 
rate for various densities. It is clearly seen from the graph that the ABP-DI model 
confined to a linear channel displays a phase transition to a polarly ordered state 
upon increase of the input power. The graph also presents a clear evidence of 
a first order phase transition, when the phase transformation is discontinuous 
(as confirmed by the standard Binder cumulant analysis). For low densities, p = 
0.04 — 0.06, we observe large jumps in the value of ip n from almost 0 in disordered 
state to 0.33 when the order is formed. We should also note that the transition 
happens earlier for more dense systems. To illustrate the phase behaviour of the 
system along the density path we plot i p n as a function of p in Fig.[5j3. The variation 
of the order parameter indicates a discontinuous phase transition on increasing the 
density. Noteworthy, the order parameter jump becomes smaller at higher densities 
and lower energy influx rates. 

The ordering in the confined system is also affected by the transverse size of the 
channel. In Fig. [6] we present the phase diagram for the ABP-DI model in channel 
confinement in coordinates p — q. We observe practically the same power law as 
in the unbounded space, q c oc p)T 0,46 HEj. In the whole range of explored densities 
and the input power, the transition is of the first order. The discontinuous nature 
of the transition is related to formation of density waves, which can be seen in the 
snapshots in Figs. [3] and [4] This issue is discussed in more detail in our paper on 
the Vicsek model [28] ■ 
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Fig. 4 Simulation snapshots (cutout of the full simulation box) of swarms for the ABP-DI 
model in a linear channel confinement at various density p (q = 0.5, T E = 0.3). (a) p = 0.04. 
(b) p = 0.1. (c) p = 0.4. Arrows indicate the direction of motion of the individuals as well as 
the magnitudes of the velocities. 

(a) (b) 




Fig. 5 Normalised order parameter for the ABP-DI system in absence of external field, Ho = 
0: (a) the iso-p curves and (b) the iso-g curves. 


3.2 Hysteresis of the mean velocity of the swarm 

To steer the swarm motion, we now apply a homogeneous oscillating field, which 
exerts a force H(t) on the particles along the channel. In simulations, we vary the 
field oscillation amplitude Hq and frequency /. The system parameters are set so 
that the swarm is orientationally ordered in absence of field. To compare different 
systems, we will further present the frequency in dimensionless form, scaled by the 
order parameter relaxation time r, which is defined by Eq. m 

Figure Q shows the measured values of the orientational order parameter as 
a function of time together with the corresponding field variation curves. At low 
frequencies, /r = 0.1, we see that the velocity of the swarm is not proportional 
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Fig. 6 Phase diagram of the ABP-DI model in absence of external field, Ho = 0, in con¬ 
finement with channel width W = 50 and in a periodic square simulation box of linear size 
L = 130. 


to the field although it changes in phase with the latter. At fr = 1 the variation 
of ip becomes sinusoidal but now exhibits a phase lag as compared to the field 
variation. At the highest frequency, fr = 10 the order parameter does not change 
the sign but rather oscillates about a fixed non-zero value \ip\ « 0.83. This value, 
however, is different from that observed in the absence of a driving field \tp\ « 0.96. 
The corresponding H — ip diagrams are presented in Fig. [8] The shape of the 
loops changes from sigmoidal at low frequency, fr = 0.1 , where the curve is also 
symmetric with respect to the change of sign of the field, to an ellipsoidal one at 
higher frequencies, fr = 0.5, 2.5 and 25. At high frequencies, as we noted before, 
the order parameter does not change the sign within the cycle, and the loops are 
not symmetric with respect to the origin. 

To illustrate the microscopic dynamic properties of the active particles, we 
plotted the instantaneous velocity histograms along with the instantaneous values 
of the field and the order parameter in Fig. [9] The motion is clearly polarised at 
fr w 0,0.25, 0.5, 0.75, and 1.0. The polarisation is strongest at fr w 0.25 and 
0.75. At the points fr = 0.1 and 0.6, where the order parameter ip(t) = 0, we 
see characteristic crater-like distributions with velocity peaks along the j/-axis. 
Therefore, the states with zero average velocity are achieved not by reduction of 
individual velocities but rather by loss of polarisation, when the particles are not 
braking to reverse the direction of motion but making an U-turn. Therefore, at 
some moments we can see them moving predominantly across the channel. 

All the loop shapes we observe here are quite familiar from the magnetic hys¬ 
teresis Wfflm- At low frequencies, where the hysteresis loop is sigmoidal, we 
measured the dynamic characteristics of the swarm, which quantify its controlla¬ 
bility. The dynamic coercivity - half-width at middle section — as can be seen in 
Fig. llOh.c.e. grows with the field frequency, field amplitude, but decreases with 
temperature. The coercivity vanishes in the static limit but grows as H c oc /°' 55 
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t/ 


Fig. 7 Variation of the order parameter for the ABP-DI swarm at three different frequencies 
of the driving field: (a) fr = 0.1. (b) /r = 1. (c) fr = 10. Other parameters: p = 0.04, 
T e = 0.3, Ho = 1. The order parameter values <p(t) are shown by the solid lines while the field 
H(t) by the dashed curves. The upper set of curves (a) corresponds to the slow field variation, 
such that the swarm has time to relax to the steady state and is always in phase with the field. 
The lower one (c) corresponds to fast field variation, such that the swarm has no chance to 
follow the field and reorient itself completely. Note that the time axis in each subplot is scaled 
by the corresponding field oscillation period. 



Fig. 8 Typical H — (/^-diagrams for the ABP-DI swarm at different frequencies of the driving 
field. Other parameters: p = 0.04, T E = 0.3, Hq = 0.4. 
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Fig. 9 Variation of the particle velocity distributions within the field oscillation cycle. System 
parameters: p = 0.04, T E = 0.3, Ho = 1, /r = 0.18. The order parameter values ip(t) are shown 
by the solid line while the field H(t) by the dashed curve. 


with the frequency. The resistance of the swarm to the action of the re-orienting 
held is related to the persistence of the particle motion and orienting action of 
the channel. In confinement, however, when the transverse size of the swarm is 
large enough, there exists a kinetic barrier for reorientation due to aligning action 
of the walls, which can prevent the reorientation and lead to long-living aligned 
states even in the presence of opposing fields. We were not able to observe static 
hysteresis at the chosen conditions, as the swarm always did reorient in the end in 
a constant held. We can envision, however, that in narrower channels and at higher 
density such that the cluster is system-spanning the system can demonstrate static 
hysteresis. In Fig. IKJb.d.f. we see that the dynamic remanence - the residual po¬ 
larisation of the swarm when the held turns zero - also grows with frequency and 
held amplitude, but decreases with temperature. In the limit of / —> 0, the dy¬ 
namic remanence simply rehects the stationary value of the mean swarm velocity 
without the held while the coercivity turns zero. At higher frequencies, however, it 
is impossible to determine these characteristics due to completely different shape 
of the loops. Both properties contribute to the integral characteristic of hysteretic 
systems, the area of the loop, which thus rehects the system’s overall dynamic con¬ 
trollability (or rather agreeability in this context) in the low frequency region. In 
a perfectly controlled system, such that the mean velocity is always in phase with 
the external held, the loop area turns zero. In contrast, a large loop area indicates 
the “amount of disagreement” between the held and the order parameter. At high 
frequencies, the held variation is too fast for the particles and their velocity hardly 
varies at all but rather oscillates around the initial value. The particles cannot 
accelerate enough. The faster the held the less they catch, so the loops ar getting 
hatter and hatter. In this region, however, it cannot be easily interpreted in terms 
of agreeability, as the particles are moving opposite to the held direction half of 
the time at any fr > 1. 
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Fig. 10 The coercivity (a),(c),(e) and the remanence (b),(d),(f) of the swarm at p = 0.04 in 
the LF region. The field amplitude is Ho = 0.4 on the temperature dependencies. The effective 
temperature is T E = 0.3 on the field dependencies, frequency is fr = 0.018 on graphs (c)-(f). 
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1 

A 


Fig. 11 The hysteresis loop area A as a function of the scaled field oscillation frequency at 
different field strengths Hq. The upper dashed curve is a fit with the Eq. GS»- Inset: Exponent 
/3 for the LF parts of the A(f) curves as a function of the field amplitude Hq. Other parameters: 
p = 0.04, T e = 0.3. 


We calculated the loop area for fixed parameters of the active particles but 
at varying field amplitude and frequencies. In Fig. 1111 we show the frequency de¬ 
pendence of the loop area at three different field magnitudes. The trends we see 
confirm our previous observations made from the shape of the loops. All the curves 
show a maximum at the reduced frequency fr m 1 and a power law decay both 
on increasing and decreasing oscillation frequency. Quite obviously, the loop area 
is greatest in the strongest field. At the smallest field amplitude, Hq = 0.1, the 
variation of the area at low frequencies is very weak. All the curves show similar 
asymptotic behaviour at /r 1: A oc 1//. At low frequencies (LF) and high fields 
the area grows proportionally to /. 

While the behaviour of the high frequency (HF) asymptotes seems to be uni¬ 
versal, the variation of the area at low frequencies is governed by power laws with 
variable exponents. In the literature on ferromagnetic materials, it is common to 
present the variation of the loop area in the form (A — Ao) oc Hq f^T~^, where Aq 
is the loop area in the static hysteresis [30]. In our model, as the system does not 
show any static hysteresis, i.e. Aq = 0, we can study the scaling of the area A as 
is. 

The exponent /? can be measured in the LF region at different system settings. 
The dependence of the exponent on the driving field amplitude Hq is shown in the 
inset in Figure fill We observe /? increase from 0.55 to 0.95 when the field grows 
from 0.1 to 10. Asymptotically, the exponent is approaching unity at Hq —> oo. 
Similarly, we can study the scaling exponent /3 at different temperatures. Figure 
1121 shows the A(/) curves at ambient temperatures from T E = 0.1 to 10. The 
exponent is growing with temperature from /? = 0.55 at T E = 0 to /3 = 0.93 at 
T e = 1. 

The sensitivity of the loop area to the field oscillation amplitude is illustrated 
in Fig. 1131 The loop area as a function of the field strength varies according to a 
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A 

Fig. 12 The hysteresis loop area A as a function of frequency / at various temperatures T E 
(p = 0.04, H 0 = 1). Inset: Exponent 3 for the LF parts of the /l(/) curves as a function of 
temperature T E . 



Fig. 13 The hysteresis loop area A variation with the field strength 11 q at various tempera¬ 
tures T e (p = 0.04, / = 0.01). Inset: Exponent a at different temperatures T^. 


power law A oc Hq with a taking values from 0.9 to 1.9 at T E varying from 0 to 
5 (see the inset). Fig. 1141 presents a similar scan along the frequency axis. Here 
we see a qualitative change of the behaviour. Firstly, all the curves show the same 
asymptotic power law A oc Hq in sufficiently high fields, Ho > 20. Secondly, in 
the I4F region, at fr > 1, this law is observed at all held strengths. However, the 
LF behaviour depends on the driving held amplitude. The exponent a grows from 
0.75 to 2 in the region fr < 1. 
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Fig. 14 The hysteresis loop area A as a function of the field strength Hq at various frequencies 
/ (p = 0.04, T e = 0.3). Inset: Exponent a as a function of scaled frequency fr. 



Fig. 15 The typical hysteresis loop area A variation with temperature T E . Inset: Exponent 
7 as a function of the field strength Hg. Other parameters: p = 0.04, / = 0.01. 


Finally, we study the influence of temperature on the loop area in the LF 
region. The simulation results are shown in Fig. 1151 and look intriguing. The loop 
area grows at low temperatures, then reaches a maximum at about T E = 0.3 and 
then shows a power law decay. The areas of the hysteresis loops decrease with the 
temperature due to the decrease of the alignment and coercivity as temperature is 
increased. The maximum does not appear at high field amplitudes. The power law 
exponent at high temperatures, as calculated from a fit with A oc {T E )~ 1 , varies 
from 7 = 1 in weak fields to 7 « 0.3 in strong fields. 
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4 Discussion 


The observations for hysteresis in the active swarm agree qualitatively with the 
corresponding results for 2D magnets [30] • First of all, we can note the special 
role of the orientational relaxation time r that determines the timescale for the 
swarm dynamics. At frequencies lower than 1/r, the active particles have enough 
time to align with the field and, therefore, the swarm follows the field direction 
“obediently”. At hight frequencies, /r » 1, the direction of motion of the particles 
does not change anymore. The action of the field leads just to a minor velocity 
oscillation about the average value: the motion is slightly slowed down in the 
opposing field but accelerated if the field is acting in the direction of motion. In 
this regime, the velocity follows the field direction only half of the time. Finally, at 
/rfs 1, the inability of the swarm to follow the direction of the held is accompanied 
by the strong variation of the velocity, so that the dynamic coercivity and dynamic 
remanence are both high, which is reflected in the large loop area. 

The variation of the loop shape and area can be understood from the follow¬ 
ing simple analysis. We consider an ensemble of non-interacting active Brownian 
particles with energy depot as described by Eqs. m and calculate an average of 
the ^-component of the acceleration over all particles: 

(^) = + c + dV? Vix + (17) 

The averaging eliminates the random term for symmetric noise. If, in addition, we 
assume a high-dissipation regime such that c S> dV 2 , we can rewrite the thrust 
term as qd/(c + dV 2 ) « qd/c — (qd 2 / c 2 )V 2 and arrive to a simple equation for 
evolution of the order parameter: 



Averaging of the second term on the right-hand side is not straightforward but 
we can assume that the term scales approximately as <p 3 [22] if the motion is 
mostly along the z-axis. The depot mode predicts that in absence of field and 
fluctuations the mean particle velocity is given by Eq. ©■ Then, in weak fields 
and for isotropic motion, the particle speed becomes independent of the direction, 

so {V 2 V lx ) » V„V 

In certain limiting cases we can derive explicit asymptotes for A(/) [32]. The 
solution can easily be found for the case when the right-hand side of Eq. m 
contains only linear term in ip and can be written as \<p + Hq sin(oit). For passive 
Brownian particles {q = 0) in a viscous medium, the second term in Eq. (1181) 
disappears while the first one is simply —”/ E ip. Moreover, the function takes the 
same form for an ABP-DI system in the strong field limit, Hq q/Vo, with 
A = -y E — dq/c playing a role of effective friction coefficient. In this case, the closed 
form for the whole A(f) curve is described by [32]: 


Mf) 


1 H 2 f 

2 (^) 2 +/ 2 ' 


(19) 


More generally, at high frequencies, fr » 1 , the time derivative of ip in Eq. ( 1181 ) be¬ 
comes very large, so ^ « Hq sinwt and a single integration gives ip oc — T cos (urt). 
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Fig. 16 The hysteresis loop area A for an ensemble of passive (q = 0) and active (q = 1) 
Brownian particles at two different field strengths: Ho = 0.2 and Ho = 1. Lines correspond 
to results obtained with non-interacting agents while symbols denote data for interacting 
particles. Other parameters: p = 0.04, T E = 0.3. 


Therefore, we find a general high frequency asymptotic result: A(f) & Hq/ 2/. In 
the strong field limit, Ho q/V6, the steady state speed of the particle is given 
by ip = (Vi x ) = Ho/^ E and the loop area scales as A oc Hq as predicted by Eq. 
(Ha. Thus, the asymptotic behaviour of A(f) at HqVq/ q 1 and /t » 1 does not 
depend on whether the particles are active or passive nor on interactions between 
them. We can summarise the limiting scaling laws as follows: 

— A oc / (d = 1) at /t -C 1 for passive Brownian particles, for non-interacting 
active particles (ksT E » D{ 7 ' E ) 3 ), or strong fields Hq » q/Vo 

— A oc f^ 1 at fr > 1 for all passive and active Brownian particles 

— A oc Hq (a = 2 ) at strong fields Ho » q/Vo 

We can clearly see these scaling laws in the simulation data presented in Figures 
1111112111411161 The low frequency law A oc f appears in Figs. 1111 and fTfil the high 
frequency one A oc / _1 - in Figs. [TT] 1121 and 1161 and the strong field asymptote 
A oc Hq is seen in Fig. [TJ] The scaling law A oc / _1 at /r > 1 seems to be valid 
at all field magnitudes, as we predicted above. 

Outside the range of these specific simple behaviours, the generic scaling form 
{A - A 0 ) oc seems to be valid. However, the exponents for the swarm 

differ from those for 2D magnets. The observed values of the exponent a (from 
0.75 to 2.0) for the swarm are higher than the numbers for the 2D Heisenberg 
model, where a = 0.40 ± 0.02 was reported [30j. Similarly, the 7 exponent at high 
temperatures, A oc taking values from 7 = 1 in weak fields to 7 « 0.3 

in strong fields, is higher than the result 2D Heisenberg model - 7 = 0.30 ± 0.02 
(30) except for the lowest value. The most obvious reason for the difference is 
the restriction on the order parameter. In the channel confinement, the order 
parameter becomes one-dimensional. In this sense, the symmetry of our problem 
is more closely resembling 2D Ising model. Indeed, in the Ising magnet, one finds 
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exponents a = 0.70, /3 = 0.36, 7 = 1.18 [ 6 ], which are close to our results at 
T e —> 0, fr —> 0, and Hq —> 0, where we have a = 0.75, 7 = 1.05. For the 
exponent /3, however, we observe values from 0.55 to 0.94, which are higher than 
/3 = 0.38 ± 0.04 for Heisenberg ferromagnet and (3 = 0.36 for Ising ferromagnet 
[6|. In our model, the exponent /3 decreases rapidly with the temperature. One 
can expect that it will reach even lower values for either lower T E or stronger 
aligning interactions. Beside the exponents, we see a qualitative analogy between 
the swarms and magnets in other properties. As was found in [30], the peak on the 
curves A(T) is observed in weak fields and is moving to lower temperatures upon 
increase of the field oscillation amplitude. In the Heisenberg magnet, however, 
the peak corresponds to a jump from non-zero loop to zero, such that below the 
critical temperature the field is unable to reorient the magnetic moment within 
the oscillation period. Our system seems to allow reorientation at any field value 
and, therefore, the jump is not observed. The reason for the area decrease at low 
temperatures is the increased resistance of a fully aligned swarm, which makes 
the reorientation within the oscillation period harder. The stronger the field, the 
shorter time it needs to reorient the swarm. That is why we do not see a peak at 
Hq = 1.0 at the chosen frequency. 

Finally, we can comment on the role of interparticle interactions and active 
propulsion in the observed dynamic hysteresis. These can be illustrated by the 
data presented in Fig. [16] '22;. For passive Brownian particles, the interactions do 
not affect the loop area: The areas obtained with or without interactions coincide 
at all frequencies. For active particles, the loop area is greater than that for passive 
ones at all frequencies as their speed is generally higher. The difference due to in¬ 
teractions is greatest at low frequencies, fr < 1. The interactions between particles 
in that region lead to the increase of the loop area, as the motion of interacting 
active swarm more is persistent than of the non-interactive one. The loop area is 
hardly affected by interactions at high frequencies, fr > 1 as the particle have no 
time to develop any collective motion. We should also note that the effect of inter¬ 
actions is most pronounced in weak fields. In strong fields, the differences between 
the interacting and non-interacting particles as well as between the passive and 
active ones are small. Thus, a study of the low frequency response and especially 
exponents a and /3, which are both always lower for the interacting systems, can 
provide information about the degree of collectivity and order in the motion of an 
active swarm. 


5 Conclusions 

We studied the dynamics of active swarms using method and ideas from condensed 
matter physics. We demonstrated that the swarms in an external field exhibit a 
dynamic hysteresis, which is qualitatively identical to that observed in magnet¬ 
ics. We measured the hysteresis loops for swarms of simulated active Brownian 
particles with dissipative interactions at various field oscillation frequencies. Our 
calculations show that the swarm reaction depends on the ratio of the orienta¬ 
tional relaxation time to the field oscillation period. At high field frequencies, the 
collective component of the behaviour becomes negligible and the swarm behaves 
as a collection of independent individuals, while at low frequencies the swarm 
can develop collective dynamics. We derived scaling exponents for several limit- 
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ing situations of the swarm dynamics: absence of interactions, weak propulsion, 
strong field, etc. All the limiting laws are confirmed by simulation results for active 
swarms. For the general case of active interacting particles, the scaling exponents 
for the hysteresis loop area are non-universal: they depend on the system’s param¬ 
eters - noise amplitude, interaction strength, and the field amplitude and do not 
coincide with exponents for 2D lattice models of ferromagnets. 
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